function c = bnet_2(dofs,V,T,d,f,varargin)
%        c = bnet_2(dofs,V,T,d,f,varargin)
% This function returns the coefficients in the vector c of the piecewise interpolation
% polynomials of the function f(x,y) over the triangulation [V,T]. That is, c is the
% coefficient vector of continuous spline function of degree d which interpolates f(x,y)
% over the domain points of each triangle in triangulation [V,T].
%
% this version assume d is nonuniform.
%
d_max = max(d);
Mat = cell(d_max,1);
for degree = 1:d_max
   Mat{degree} = vdm21(degree);
end
m = 0; nt = size(T,1);
for k=1:nt
    degree = d(k);
    md = (degree+1)*(degree+2)/2;
    m = m + md^2;
end
Indx1 = zeros(m,1);
Indx2 = zeros(m,1);
S = zeros(m,1);
dim = max(max(dofs));
b = zeros(dim,1);
pos_start = 1;
for tri = 1:nt
    [I,J,K]=indices(d(tri));
    pts = (I*V(T(tri,1),:)+J*V(T(tri,2),:)+K*V(T(tri,3),:))/d(tri);
    b_loc = feval(f,pts(:,1),pts(:,2),varargin{:});  
    loc_dof = dofs(tri,1):dofs(tri,2);
    mat = Mat{d(tri)};
    [i,j,s] = find(mat);
    L = length(i);
    Indx1(pos_start:(pos_start + L-1)) = loc_dof(i);
    Indx2(pos_start:(pos_start + L-1)) = loc_dof(j);
    S(pos_start:(pos_start + L-1)) = s;
    pos_start = pos_start + L;        
    b(loc_dof) = b(loc_dof) + b_loc;
end;
A = sparse(Indx1(1:(pos_start-1)),Indx2(1:(pos_start-1)),S(1:(pos_start-1)),dim,dim);
c = A\b;
